auto_annot_Haber2017_with_Smillie2019_dblabel_level3

import besca as bc
import scanpy as sc
import pandas as pd
import pkg_resources
import os

Specify folders where .h5ad files are found and their names.

The datasets that are already annotated and should be used for training. If you only use one dataset please use list of one.

adata_trains = [bc.datasets.Smillie2019_processed()]

The dataset of interest that should be annotated.

adata_pred = bc.datasets.Haber2017_processed()
adata_orig = bc.datasets.Haber2017_processed()

Define level of dblabel reference annotation

level = 3

Give your analysis a name.

analysis_name = 'auto_annot_Haber2017_with_Smillie2019_dblabel_l' + str(level) 

Now specify parameters

Specify column name of celltype annotation you want to train on.

celltype_train ='dblabel_l' + str(level)
celltype_test = 'dblabel_l' + str(level)

Choose a method:

• linear: Support Vector Machine with Linear Kernel
• sgd: Support Vector Machine with Linear Kernel using Stochastic Gradient Descent
• rbf: Support Vector Machine with radial basis function kernel. Very time intensive, use only on small datasets.
• logistic_regression: Standard logistic classifier iwth multinomial loss.
• logistic_regression_ovr: Logistic Regression with one versus rest classification.
• logistic_regression_elastic: Logistic Regression with elastic loss, cross validates among multiple l1 ratios.

method = 'logistic_regression_elastic' # 'logistic_regression'

Specify merge method if using multiple training datasets. Needs to be either scanorama or naive.

merge = 'scanorama'

Decide if you want to use the raw format or highly variable genes. Raw increases computational time and does not necessarily improve predictions.

use_raw = False

You can choose to only consider a subset of genes from a signature set.

genes_to_use = 'all'

Translate cell type annoation to lower dblabel level

new_cnames = bc.tl.sig.obtain_new_label(
    nomenclature_file=pkg_resources.resource_filename('besca', 'datasets/nomenclature/CellTypes_v1.tsv'), 
    cnames=list(adata_trains[0].obs['dblabel'].cat.categories), 
    reference_label='dblabel', 
    new_label='dblabel', 
    new_level=level)

new_cnames

	new_label
CD1c-positive myeloid dendritic cell	myeloid dendritic cell
CD4-positive, alpha-beta memory T cell	T cell
CD8-positive, alpha-beta T cell	T cell
CD8-positive, alpha-beta cytokine secreting effector T cell	T cell
CD141-positive myeloid dendritic cell	myeloid dendritic cell
HEV endothelial cell	HEV endothelial cell
activated CD4-positive, alpha-beta T cell	T cell
brush cell	brush cell
endothelial cell	endothelial cell
enterocyte	intestinal epithelial cell
enterocyte progenitor	intestinal epithelial cell
enteroendocrine cell	enteroendocrine cell
exhausted-like CD4-positive, alpha-beta T cell	T cell
fibroblast	fibroblast
follicular B cell	B cell
germinal center B cell	B cell
glial cell	glial cell
goblet cell	goblet cell
immature enterocyte	intestinal epithelial cell
immature goblet cell	goblet cell
inflammatory fibroblast	inflammatory fibroblast
inflammatory monocyte	monocyte
innate lymphoid cell	innate lymphoid cell
macrophage	macrophage
mast cell	mast cell
microfold cell	microfold cell
microvascular endothelial cell	blood vessel endothelial cell
myofibroblast cell	myofibroblast cell
natural killer cell	innate lymphoid cell
pericyte cell	pericyte cell
plasma cell	plasma cell
proliferating B cell	B cell
proliferating T cell	proliferating T cell
proliferating monocyte	proliferating monocyte
proliferating transit amplifying cell	intestinal epithelial cell
regulatory T cell	T cell
stem cell	stem cell
transit amplifying cell	intestinal epithelial cell

adata_trains[0].obs['dblabel_l' + str(level)] = bc.tl.sig.add_anno(adata_trains[0], new_cnames, 'new_label', 'dblabel')

new_cnames = bc.tl.sig.obtain_new_label(
    nomenclature_file=pkg_resources.resource_filename('besca', 'datasets/nomenclature/CellTypes_v1.tsv'), 
    cnames=list(adata_pred.obs['dblabel'].cat.categories), 
    reference_label='dblabel', 
    new_label='dblabel', 
    new_level=level)

new_cnames

	new_label
goblet cell	goblet cell
proliferating epithelial fate stem cell	somatic stem cell
enterocyte	intestinal epithelial cell
epithelial fate stem cell	somatic stem cell
immature goblet cell	goblet cell
brush cell	brush cell
proliferating transit amplifying cell	intestinal epithelial cell
paneth cell	intestinal epithelial cell
transit amplifying cell	intestinal epithelial cell
proliferating enterocyte progenitor	intestinal epithelial cell
immature enterocyte	intestinal epithelial cell
enteroendocrine cell	enteroendocrine cell
enterocyte progenitor	intestinal epithelial cell

adata_pred.obs['dblabel_l' + str(level)] = bc.tl.sig.add_anno(adata_pred, new_cnames, 'new_label', 'dblabel')
adata_orig.obs['dblabel_l' + str(level)] = adata_pred.obs['dblabel_l' + str(level)] 

Prepare all training and the testing set.

# Select epithelial subset from Smillie2019 dataset
epithelial_subset = bc.subset_adata(adata_trains[0], adata_trains[0].obs.celltype_highlevel == 'Epi', raw=False)
adata_trains[0] = epithelial_subset

# Convert mouse symbols (MGI) to human symbols (HGNC)
mousehuman_file = pkg_resources.resource_filename('besca', 'datasets/homologs/MGItoHGNC.csv')
mousehuman=pd.read_csv(mousehuman_file,sep='\t',header='infer', encoding="unicode_escape")
mousehuman.index=mousehuman['MGI']
conversion=pd.Series(data=mousehuman['HGNC'], index=mousehuman.index)

# Convert mouse symbols (MGI) to human symbols (HGNC)
adata_orig.var.rename(columns={'SYMBOL':'MGI'}, inplace=True)
adata_orig.var['SYMBOL'] = adata_orig.var['MGI'].map(lambda x: conversion.get(x, default='') if type(conversion.get(x, default='')) == str else conversion.get(x, default=None).values[0])
adata_orig.var.index = adata_orig.var.SYMBOL
adata_orig.var_names_make_unique()
adata_pred = adata_orig.copy()

This function merges training datasets, removes unwanted genes, and if scanorama is used corrects for datasets.

adata_train, adata_pred = bc.tl.auto_annot.merge_data(adata_trains, adata_pred, genes_to_use = genes_to_use, merge = merge)

merging with scanorama
using scanorama rn
Found 278 genes among all datasets
[[0.         0.75137665]
 [0.         0.        ]]
Processing datasets (0, 1)
integrating training set
calculating intersection

Train the classifier.

The returned scaler is fitted on the training dataset (to zero mean and scaled to unit variance).

classifier, scaler = bc.tl.auto_annot.fit(adata_train, method, celltype_train)

[Parallel(n_jobs=10)]: Using backend ThreadingBackend with 10 concurrent workers.

max_iter reached after 63 seconds
max_iter reached after 63 seconds
max_iter reached after 63 seconds
max_iter reached after 63 seconds
max_iter reached after 66 seconds
max_iter reached after 67 seconds
max_iter reached after 68 seconds
max_iter reached after 68 seconds
max_iter reached after 68 seconds
max_iter reached after 65 seconds
max_iter reached after 65 seconds
max_iter reached after 65 seconds
max_iter reached after 65 seconds
max_iter reached after 65 seconds
max_iter reached after 65 seconds
max_iter reached after 64 seconds
max_iter reached after 65 seconds
max_iter reached after 65 seconds
max_iter reached after 67 seconds
max_iter reached after 68 seconds
max_iter reached after 68 seconds
max_iter reached after 67 seconds
max_iter reached after 68 seconds
max_iter reached after 67 seconds
max_iter reached after 67 seconds
max_iter reached after 66 seconds
max_iter reached after 67 seconds
max_iter reached after 65 seconds
max_iter reached after 65 seconds
max_iter reached after 66 seconds

[Parallel(n_jobs=10)]: Done   3 out of   3 | elapsed: 11.0min finished

Prediction

Use fitted model to predict celltypes in adata_pred. Prediction will be added in a new column called 'auto_annot'. Paths are needed as adata_pred will revert to its original state (all genes, no additional corrections). The threshold should be set to 0 or left out for SVM. For logisitic regression the threshold can be set.

adata_predicted = bc.tl.auto_annot.adata_predict(classifier = classifier, scaler = scaler, adata_pred = adata_pred, adata_orig = adata_orig, threshold = 0)

Write out metrics to a report file, create confusion matrices and comparative umap plots

adata_pred.obs

	CELL	CONDITION	sample_type	donor	region_x	sample	percent_mito	n_counts	n_genes	batch	...	celltype0	celltype1	celltype2	celltype3	dblabel	barcode	region_y	cell_label	_merge	dblabel_l3
0	haber_intestine_donor_M1_Duo.AAACATACAGCGGA	healthy	mouse_small_intestine_epithelial	M1	Duo	Duo_M1	0.001410	12768.0	1227	Duo	...	epithelial cell	paneth cell	paneth cell	paneth cell	paneth cell	AAACATACAGCGGA	Jejunum	Paneth	both	intestinal epithelial cell
1	haber_intestine_donor_M1_Duo.AAACATACCTTACT	healthy	mouse_small_intestine_epithelial	M1	Duo	Duo_M1	0.010779	6583.0	2156	Duo	...	epithelial cell	enterocyte	enterocyte	enterocyte	enterocyte	AAACATACCTTACT	Jejunum	Enterocyte	both	intestinal epithelial cell
2	haber_intestine_donor_M1_Duo.AAACCGTGCAGTCA	healthy	mouse_small_intestine_epithelial	M1	Duo	Duo_M1	0.022508	2799.0	1362	Duo	...	epithelial cell	epithelial fate stem cell	proliferating epithelial fate stem cell	proliferating epithelial fate stem cell	proliferating epithelial fate stem cell	AAACCGTGCAGTCA	Jejunum	TA	both	somatic stem cell
3	haber_intestine_donor_M1_Duo.AAACGCTGCAGTCA	healthy	mouse_small_intestine_epithelial	M1	Duo	Duo_M1	0.015041	6048.0	2287	Duo	...	epithelial cell	transit amplifying cell	transit amplifying cell	transit amplifying cell	transit amplifying cell	AAACGCTGCAGTCA	Jejunum	TA	both	intestinal epithelial cell
4	haber_intestine_donor_M1_Duo.AAACGCTGCGTGAT	healthy	mouse_small_intestine_epithelial	M1	Duo	Duo_M1	0.023022	2780.0	1320	Duo	...	epithelial cell	transit amplifying cell	transit amplifying cell	transit amplifying cell	transit amplifying cell	AAACGCTGCGTGAT	Jejunum	Stem	both	intestinal epithelial cell
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
10891	haber_intestine_donor_M2_Il.TTTCAGTGACCAGT	healthy	mouse_small_intestine_epithelial	M2	Il	Il_M2	0.008591	5468.0	2037	Il	...	epithelial cell	transit amplifying cell	transit amplifying cell	transit amplifying cell	transit amplifying cell	TTTCAGTGACCAGT	Ileum	EP	both	intestinal epithelial cell
10892	haber_intestine_donor_M2_Il.TTTCGAACAGAACA	healthy	mouse_small_intestine_epithelial	M2	Il	Il_M2	0.007760	10174.0	2884	Il	...	epithelial cell	immature enterocyte	immature enterocyte	immature enterocyte	immature enterocyte	TTTCGAACAGAACA	Ileum	Enterocyte	both	intestinal epithelial cell
10893	haber_intestine_donor_M2_Il.TTTCTACTGCTCCT	healthy	mouse_small_intestine_epithelial	M2	Il	Il_M2	0.006121	9307.0	2856	Il	...	epithelial cell	immature goblet cell	immature goblet cell	immature goblet cell	immature goblet cell	TTTCTACTGCTCCT	Ileum	Goblet	both	goblet cell
10894	haber_intestine_donor_M2_Il.TTTGACTGCGCCTT	healthy	mouse_small_intestine_epithelial	M2	Il	Il_M2	0.007189	4029.0	1277	Il	...	epithelial cell	goblet cell	goblet cell	goblet cell	goblet cell	TTTGACTGCGCCTT	Ileum	Goblet	both	goblet cell
10895	haber_intestine_donor_M2_Il.TTTGCATGGAGGAC	healthy	mouse_small_intestine_epithelial	M2	Il	Il_M2	0.009316	8907.0	2277	Il	...	epithelial cell	enterocyte	enterocyte	enterocyte	enterocyte	TTTGCATGGAGGAC	Ileum	Enterocyte	both	intestinal epithelial cell

10896 rows × 21 columns

%matplotlib inline

bc.tl.report(
    adata_pred=adata_predicted,
    celltype=celltype_test, 
    method=method, 
    analysis_name=analysis_name,
    train_datasets = adata_trains,
    test_dataset = adata_orig,
    merge = merge,
    name_prediction='auto_annot', 
    name_report='auto_annot', 
    use_raw=use_raw,
    remove_nonshared=True, 
    clustering='leiden', 
    asymmetric_matrix=True, 
    delimiter='\t', 
    verbose=True
)

acc: 0.57
f1: 0.52

... storing 'auto_annot' as categorical
... storing 'SYMBOL' as categorical

ami: 0.35
ari: 0.26
silhouette dblabel_l3: 0.24
silhouette auto_annot: 0.05
pair confusion matrix:
          0         1
0  42139704  33801324
1  11555152  31215740
WARNING: saving figure to file figures/umap.ondata_auto_annot_Haber2017_with_Smillie2019_dblabel_l3.png

WARNING: saving figure to file figures/umap.auto_annot_Haber2017_with_Smillie2019_dblabel_l3_dblabel_l3.png

WARNING: saving figure to file figures/umap.auto_annot_Haber2017_with_Smillie2019_dblabel_l3_auto_annot.png

WARNING: saving figure to file figures/umap.auto_annot_Haber2017_with_Smillie2019_dblabel_l3_leiden.png

sc.pl.umap(adata_predicted, color=[celltype_test, 'auto_annot'])
sc.pl.umap(adata_predicted, color=[celltype_test, 'auto_annot'], legend_loc='on data', legend_fontsize=8)

sc.pl.umap(adata_predicted, color=[celltype_test, 'auto_annot'], legend_fontsize=7, wspace = 1.4, save = '.svg')
sc.pl.umap(adata_predicted, color=[celltype_test, 'auto_annot'], legend_loc='on data', legend_fontsize=7, wspace = 1.4, save = '.ondata.svg')

WARNING: saving figure to file figures/umap.svg

WARNING: saving figure to file figures/umap.ondata.svg

adata_train

View of AnnData object with n_obs × n_vars = 46102 × 278
    obs: 'CELL', 'Cluster', 'Health', 'Location', 'Subject', 'celltype_highlevel', 'nGene', 'nUMI', 'original_name', 'percent_mito', 'n_counts', 'n_genes', 'batch', 'leiden', 'dblabel', 'celltype', 'cluster_celltype', 'Type', 'dblabel_l3'
    var: 'SYMBOL', 'ENSEMBL-0', 'n_cells-0', 'total_counts-0', 'frac_reads-0', 'ENSEMBL-1', 'n_cells-1', 'total_counts-1', 'frac_reads-1', 'MGI-1', 'highly_variable-1', 'means-1', 'dispersions-1', 'dispersions_norm-1', 'mean-1', 'std-1'
    uns: 'Cluster_colors', 'Location_colors', 'Type_colors', 'celltype_highlevel_colors', 'leiden', 'leiden_colors', 'umap'
    obsm: 'X_pca', 'X_umap', 'X_scanorama'

adata_predicted_wo_unknown = adata_predicted.copy()
adata_predicted_wo_unknown = bc.subset_adata(adata_predicted_wo_unknown, adata_predicted_wo_unknown.obs.auto_annot != 'unknown', raw=False)
bc.pl.riverplot_2categories(adata_predicted_wo_unknown, [celltype_test, 'auto_annot'])

# Compare to random assignment
import random
random.seed(1)
adata_predicted.obs['random_labeling'] = list(adata_predicted.obs[celltype_test].sample(frac=1))

bc.tl.report(
    adata_pred=adata_predicted,
    celltype=celltype_test, 
    method="compare_to_random_" + method, 
    analysis_name=analysis_name,
    train_datasets = adata_trains,
    test_dataset = adata_orig,
    merge = merge,
    name_prediction="random_labeling",
    name_report="compare_to_random_auto_annot", 
    use_raw=use_raw,
    remove_nonshared=False, 
    clustering='leiden', 
    asymmetric_matrix=True, 
    delimiter='\t', 
    verbose=True)

acc: 0.36
f1: 0.36

... storing 'random_labeling' as categorical

ami: 0.0
ari: 0.0
silhouette dblabel_l3: 0.24
silhouette random_labeling: -0.02
pair confusion matrix:
          0         1
0  48662900  27278128
1  27278128  15492764
WARNING: saving figure to file figures/umap.ondata_auto_annot_Haber2017_with_Smillie2019_dblabel_l3.png

WARNING: saving figure to file figures/umap.auto_annot_Haber2017_with_Smillie2019_dblabel_l3_dblabel_l3.png

WARNING: saving figure to file figures/umap.auto_annot_Haber2017_with_Smillie2019_dblabel_l3_random_labeling.png

WARNING: saving figure to file figures/umap.auto_annot_Haber2017_with_Smillie2019_dblabel_l3_leiden.png

from sinfo import sinfo
sinfo()

-----
anndata             0.7.5
besca               2.4+57.g5ad53b2
pandas              1.2.2
pkg_resources       NA
plotly              4.14.3
scanpy              1.6.1
sinfo               0.3.1
sklearn             0.24.1
-----
IPython             7.20.0
jupyter_client      6.1.11
jupyter_core        4.7.1
notebook            6.2.0
-----
Python 3.7.9 | packaged by conda-forge | (default, Dec  9 2020, 21:08:20) [GCC 9.3.0]
Linux-3.10.0-693.11.6.el7.x86_64-x86_64-with-centos-7.4.1708-Core
24 logical CPU cores, x86_64
-----
Session information updated at 2021-07-18 09:08

Convert to html

%%javascript

IPython.notebook.kernel.execute('nb_name = "' + IPython.notebook.notebook_name + '"')

nb_name = os.path.join(os.getcwd(), nb_name)

! jupyter nbconvert --to html {nb_name}

[NbConvertApp] Converting notebook /pstore/data/bioinfo/users/hatjek/devel/besca_publication_results/intestine/auto_annot/auto_annot_Haber2017_with_Smillie2019_dblabel_l3.ipynb to html
[NbConvertApp] Writing 6902248 bytes to /pstore/data/bioinfo/users/hatjek/devel/besca_publication_results/intestine/auto_annot/auto_annot_Haber2017_with_Smillie2019_dblabel_l3.html